{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install pandas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from scipy import stats\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a data scientist you'll spend a lot of time answering questions with data. I currently work as a data scientist in the healthcare industry, providing metrics and building models for hospitals and healthcare related organizations. In my practice most of my time is spent doing two things:\n",
    "* Translating qualitative business questions into rigorous solutions that can be generated with data\n",
    "* Implementing these solutions programmatically\n",
    "\n",
    "I'm going to walk you through two questions that I've actually been asked on the job:\n",
    "1. Should my hospital focus more on improving its mortality rate?\n",
    "2. Which pharmacists are handing out too many opioids?\n",
    "\n",
    "What they both have in common is that there's a right way and a wrong way to do them. Further, it's very easy to answer these questions the wrong way and have your mistakes go unnoticed. As you'll see, the difference between great and average answers to these questions is a having a little bit of a mathematical background."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem Statement 1:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Centers for Medicare & Medicaid Services (CMS) is responsible for rating hospitals based off of quality and performance metrics. \n",
    "\n",
    "They provide various metrics to estimate hospital performance and one overarching metric known as the `quality_rating`. Another less important rating is the `mortality_rate`. \n",
    "\n",
    "Suppose you're tasked with investigating the relationship between a hospital's `quality_rating` and its `mortality_rate` to answer the first question I mentioned above. How much does a hospital's mortality rate affect its overall `quality_rating`? Assume the following about your data:\n",
    "\n",
    "* `quality_rating`: A numerical value in the range (1, 2, 3, 4, 5)\n",
    "* `mortality_rate`: One of the following values:\n",
    "    * *above_average*\n",
    "    * *average*\n",
    "    * *below_average*\n",
    "    \n",
    "And assume you're given the data below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of Hospitals: 100\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mortality_rate</th>\n",
       "      <th>quality_rating</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>average</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>below_average</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>above_average</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>average</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>average</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>below_average</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>below_average</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>below_average</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>average</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>average</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  mortality_rate  quality_rating\n",
       "0        average               1\n",
       "1  below_average               4\n",
       "2  above_average               1\n",
       "3        average               4\n",
       "4        average               4\n",
       "5  below_average               5\n",
       "6  below_average               4\n",
       "7  below_average               4\n",
       "8        average               3\n",
       "9        average               5"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ratings_df = pd.read_csv('hospital_ratings.csv').drop(columns='Unnamed: 0')\n",
    "\n",
    "print('Number of Hospitals: {}'.format(len(ratings_df.index)))\n",
    "\n",
    "ratings_df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, each row represents a hospital. Each hospital has a mortality rate and quality score. The premise is to *compare* `mortality_rate` and `quality_rating` so a key statistical tool should be jumping to your head. **Correlation**. \n",
    "\n",
    "This is the first step of reducing our *qualitative* question into *quantitative* terms. I urge you to take a stab at the problem at this point — or at least flesh out a complete plan to answering the question: \n",
    "\n",
    "***What's the correlation between these two metrics?***\n",
    "\n",
    "Did you arrive at a solution? Lets walk through this together."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "#### Encoding mortality rate\n",
    "\n",
    "There are arguably two steps remaining in this problem. First, we should note that our `mortality_rate` column contains ordinal data. We can only perform correlation on numerical data so we'll have to do something about this. \n",
    "\n",
    "Let's encode our data in some way that preserves its ordinal relationship. In other words, we want to our encoding for ***below average*** `mortality_rate` to be numerically superior to our encoding for ***above average*** `mortality_rate` (Note that a hospital having a **lower** mortality rate is better).\n",
    "\n",
    "How did you tackle encoding your data? There's an infinite amount of options:\n",
    "\n",
    "For instance:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "    'below average' --> 5\n",
    "    'average' --> 3\n",
    "    'above average' --> 1\n",
    "    \n",
    "or:\n",
    "\n",
    "    'below average' --> 3\n",
    "    'average' --> 2\n",
    "    'above average' --> 1\n",
    "    \n",
    "or maybe:\n",
    "\n",
    "    'below average' --> 1\n",
    "    'average' --> 0\n",
    "    'above average' --> -1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you made the decision to encode your data, what option did you settle on? Which one is correct? The correct answer is *none* of them — or *all* of them I guess...\n",
    "\n",
    "I'll explain why in a bit. But for now, let’s just continue with mapping to last option: mapping to the set: (-1, 0, 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mortality_rate</th>\n",
       "      <th>quality_rating</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   mortality_rate  quality_rating\n",
       "0               0               1\n",
       "1               1               4\n",
       "2              -1               1\n",
       "3               0               4\n",
       "4               0               4"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MORTALITY_MAP = {\n",
    "    'below_average': 1,\n",
    "    'average': 0,\n",
    "    'above_average': -1\n",
    "}\n",
    "\n",
    "ratings_df['mortality_rate'] = ratings_df['mortality_rate'].map(MORTALITY_MAP) # apply the map\n",
    "ratings_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last step to this problem is (deceptively) the simplest. Now we just need to correlate our two columns right?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mortality_rate</th>\n",
       "      <th>quality_rating</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>mortality_rate</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.473712</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>quality_rating</th>\n",
       "      <td>0.473712</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                mortality_rate  quality_rating\n",
       "mortality_rate        1.000000        0.473712\n",
       "quality_rating        0.473712        1.000000"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ratings_df.corr()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not quite. We've made a huge assumption about our data. We've assumed that our data was anything *more* than ordinal.\n",
    "\n",
    "#### Non-parametric correlation\n",
    "\n",
    "As a data scientist, you should be familiar with the concept of [parametric](https://en.wikipedia.org/wiki/Parametric_statistics) and [non-parametric](https://en.wikipedia.org/wiki/Nonparametric_statistics) statistics. Most statistical correlations will default to a parametric correlation method known as `pearson`. In our case, we need to use a non-parametric correlation known as `spearman`. \n",
    "\n",
    "This non-parametric method assumes that our data is nothing more than ordinal. In other words, no matter what mapping we make, the spearman correlation will return the same result — so long as the order is preserved after mapping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/3bcYLj11uME\" frameborder=\"0\" allow=\"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%HTML\n",
    "<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/3bcYLj11uME\" frameborder=\"0\" allow=\"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we specify that we need to use the spearman correlation method, then we notice a drastic change in results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mortality_rate</th>\n",
       "      <th>quality_rating</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>mortality_rate</th>\n",
       "      <td>1.000000</td>\n",
       "      <td>0.379098</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>quality_rating</th>\n",
       "      <td>0.379098</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                mortality_rate  quality_rating\n",
       "mortality_rate        1.000000        0.379098\n",
       "quality_rating        0.379098        1.000000"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ratings_df.corr(method='spearman')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've just dropped the correlation by about **10%**. That's a drastic difference! \n",
    "\n",
    "Had you reported this correlation using pearson as opposed to spearman, you might have seriously mislead a customer or coworker. \n",
    "\n",
    "This goes to show that having a **strong statistical background is essential for the data science role**. Such a deceptively simple problem had a key mathematical step that's often overlooked."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem Statement 2:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have millions of pharmaceutical records for the past 5 years in Massachusetts, a state well known for its opioid abuse problem. For this reason, it would be very useful if you could identify pharmacists who are prescribing too many opioids.\n",
    "\n",
    "<figure>\n",
    "    <img src=\"https://d14rmgtrwzf5a.cloudfront.net/sites/default/files/ratemass.jpg\"/>\n",
    "    <figcaption>*Source: [National Institute on Drug Abuse](https://www.drugabuse.gov/drugs-abuse/opioids/opioid-summaries-by-state/massachusetts-opioid-summary)*</figcaption>\n",
    "</figure>\n",
    "\n",
    "You're tasked with the following: \n",
    "\n",
    "***Generate a list of pharmacists who are overprescribing opioids.***\n",
    "\n",
    "Seems simple enough. Suppose you have the following dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>prescriber_id</th>\n",
       "      <th>num_opioid_prescriptions</th>\n",
       "      <th>num_prescriptions</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>BDRSDDB2&amp;RA</td>\n",
       "      <td>4</td>\n",
       "      <td>7</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3$N3WG4&amp;BBO</td>\n",
       "      <td>25</td>\n",
       "      <td>28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>S2NGGOGBRB4</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>S2NGGOGBRB4</td>\n",
       "      <td>18</td>\n",
       "      <td>22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NBORABOESOI</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>4O4EAGWROGB</td>\n",
       "      <td>20</td>\n",
       "      <td>30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>WPOARS4GW9I</td>\n",
       "      <td>11</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>PDS$GG9GOGA</td>\n",
       "      <td>29</td>\n",
       "      <td>37</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>BBP34DOGADB</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>S$BO&amp;OGO99O</td>\n",
       "      <td>16</td>\n",
       "      <td>19</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  prescriber_id  num_opioid_prescriptions  num_prescriptions\n",
       "0   BDRSDDB2&RA                         4                  7\n",
       "1   3$N3WG4&BBO                        25                 28\n",
       "2   S2NGGOGBRB4                         1                  1\n",
       "3   S2NGGOGBRB4                        18                 22\n",
       "4   NBORABOESOI                         1                  1\n",
       "5   4O4EAGWROGB                        20                 30\n",
       "6   WPOARS4GW9I                        11                 12\n",
       "7   PDS$GG9GOGA                        29                 37\n",
       "8   BBP34DOGADB                         1                  1\n",
       "9   S$BO&OGO99O                        16                 19"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prescribers_df = pd.read_csv('prescriptions.csv').drop(columns='Unnamed: 0')\n",
    "prescribers_df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each row represents a days worth of prescriptions\n",
    "\n",
    "And we have the following columns:\n",
    "* `prescriber_id`: A random code that identifies pharmacists\n",
    "* `num_opioid_prescriptions`: The number of opioid prescriptions prescribed on a given day\n",
    "* `num_prescriptions`: The number of *total* prescriptions prescribed on a given day\n",
    "\n",
    "### Grouping prescribers\n",
    "The first thing we should note is that `prescriber_id` does not contain unique values as prescribers may have prescribed drugs over multiple days. Since this metric is on a prescriber level — not a daily level — we should fix that with a pandas `groupby`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>num_opioid_prescriptions</th>\n",
       "      <th>num_prescriptions</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>prescriber_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>$39AGIIOGO2</th>\n",
       "      <td>15</td>\n",
       "      <td>22</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$DG2GOGDD9D</th>\n",
       "      <td>14</td>\n",
       "      <td>18</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$OBDNSADD96</th>\n",
       "      <td>37</td>\n",
       "      <td>57</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>&amp;$4OO2EW$OB</th>\n",
       "      <td>16</td>\n",
       "      <td>20</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>&amp;64O2OASRNA</th>\n",
       "      <td>22</td>\n",
       "      <td>23</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>&amp;AGO$DEOPRR</th>\n",
       "      <td>47</td>\n",
       "      <td>58</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>&amp;OS2DN4DADO</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>&amp;SOEB4RODOG</th>\n",
       "      <td>39</td>\n",
       "      <td>51</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2BDPWROOBOR</th>\n",
       "      <td>13</td>\n",
       "      <td>17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3$N3WG4&amp;BBO</th>\n",
       "      <td>32</td>\n",
       "      <td>39</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               num_opioid_prescriptions  num_prescriptions\n",
       "prescriber_id                                             \n",
       "$39AGIIOGO2                          15                 22\n",
       "$DG2GOGDD9D                          14                 18\n",
       "$OBDNSADD96                          37                 57\n",
       "&$4OO2EW$OB                          16                 20\n",
       "&64O2OASRNA                          22                 23\n",
       "&AGO$DEOPRR                          47                 58\n",
       "&OS2DN4DADO                           1                  1\n",
       "&SOEB4RODOG                          39                 51\n",
       "2BDPWROOBOR                          13                 17\n",
       "3$N3WG4&BBO                          32                 39"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prescribers_df = prescribers_df.groupby('prescriber_id').agg('sum')\n",
    "prescribers_df.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Much better. Now the task at hand is to rank prescribers on the basis of how many opioids they prescribe. But prescribers prescribe different amounts of drugs. This hints that we should now take each prescriber's opioid prescription ratio. Let's do just that.\n",
    "\n",
    "### Opioid prescription ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "prescribers_df['opioid_prescription_ratio'] = (prescribers_df['num_opioid_prescriptions'] / \n",
    "                                               prescribers_df['num_prescriptions'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pandas made that extremely easy. So we're done right? Just sort the prescribers by `opioid_ratio` and tie this task off?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>num_opioid_prescriptions</th>\n",
       "      <th>num_prescriptions</th>\n",
       "      <th>opioid_prescription_ratio</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>prescriber_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>DE$W$6WBN9W</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3NG4ONN46GO</th>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6W349S&amp;GD&amp;N</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4WA6NBOO&amp;DB</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>NBORABOESOI</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4R9RDN2G6D9</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>OAIGB9OGSE$</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PGWSD49S6OW</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>BO&amp;2BGDNGO&amp;</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4BBA2OGP9G6</th>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               num_opioid_prescriptions  num_prescriptions  \\\n",
       "prescriber_id                                                \n",
       "DE$W$6WBN9W                           1                  1   \n",
       "3NG4ONN46GO                           2                  2   \n",
       "6W349S&GD&N                           1                  1   \n",
       "4WA6NBOO&DB                           1                  1   \n",
       "NBORABOESOI                           1                  1   \n",
       "4R9RDN2G6D9                           1                  1   \n",
       "OAIGB9OGSE$                           1                  1   \n",
       "PGWSD49S6OW                           1                  1   \n",
       "BO&2BGDNGO&                           1                  1   \n",
       "4BBA2OGP9G6                           1                  1   \n",
       "\n",
       "               opioid_prescription_ratio  \n",
       "prescriber_id                             \n",
       "DE$W$6WBN9W                          1.0  \n",
       "3NG4ONN46GO                          1.0  \n",
       "6W349S&GD&N                          1.0  \n",
       "4WA6NBOO&DB                          1.0  \n",
       "NBORABOESOI                          1.0  \n",
       "4R9RDN2G6D9                          1.0  \n",
       "OAIGB9OGSE$                          1.0  \n",
       "PGWSD49S6OW                          1.0  \n",
       "BO&2BGDNGO&                          1.0  \n",
       "4BBA2OGP9G6                          1.0  "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prescribers_df.sort_values('opioid_prescription_ratio', ascending=False).head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not quite. This data should raise some red flags for you.\n",
    "\n",
    "Does it really make sense to rank a pharmacist as an abuser for prescribing 1 opioid out of 1 of their total prescriptions? Try to put into words exactly what's wrong with this situation. Why do we not want to report these prescribers as the worst offenders? Hopefully the answer you arrive at is something along the lines of: \n",
    "\n",
    "***Because we don't have enough information to judge them.***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imagine you were sitting in a pharmacy watching a pharmacist (let's call them \"Bill\") perform their job. Suppose you want to report Bill if he prescribes too many opioids. \n",
    "\n",
    "A customer walks in and Bill prescribes them *hydrocodone* (a common opioid). It'd be silly to report Bill immediately for this. You'd want to wait and see more of Bill's actions before you judged him because you're not confident in your findings yet. \n",
    "\n",
    "**Confidence** is a key word here as we'll soon find out."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Building confidence\n",
    "\n",
    "So how do we fix this issue? \n",
    "\n",
    "Statistically speaking, prescribing an opioid or not can be though of as a [Bernouli parameter](https://en.wikipedia.org/wiki/Bernoulli_distribution) (a fancy term for something that's binary in value — either true or false). And given the number of observations of a bernoulli parameter, we want to predict this parameter's true value. \n",
    "\n",
    "By \"true value\" I mean the actual value that prescriber's `opioid_prescription_ratio` would converge to if we had enough observations. The \"true value\" in our example with Bill would be equivalent to Bill's `opioid_prescription_ratio` if we were able to watch and record his actions for a very long time — like a year."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you've taken a statistics course you're probably familiar with the concept of confidence intervals. If not, in simplest terms, a confidence interval is just a measure of mathematical confidence attached to a range that you believe an unknown value exists within. \n",
    "\n",
    "<figure>\n",
    "    <img src=\"https://cdn-images-1.medium.com/max/2000/1*qSCzTfliGMCcPfIQcGIAJw.jpeg\" width=\"700px\" />\n",
    "    <figcaption>*Source: [Who needs backup dancers when you can have confidence intervals?](https://medium.com/design-ibm/who-needs-backup-dancers-when-you-can-have-confidence-intervals-485f9464c06f)*</figcaption>\n",
    "</figure>\n",
    "\n",
    "For example, I could be 95% confident that tomorrow's temperature will fall between 40 Fahrenheit and 70 Fahrenheit. \n",
    "\n",
    "In 1927 mathematician Edwin Wilson applied the concept of confidence intervals to Bernoulli parameters. This means we can guess the *true* value of a pharmacists `opioid_prescription_ratio` given the amount of data we have on them! \n",
    "\n",
    "Here's the formula:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wilson Confidence Interval Lower Bound: \n",
    "\n",
    "$$\\frac{\\hat{p}+\\frac{z^2}{2n}-z\\sqrt{[\\hat{p}(1-\\hat{p})+\\frac{z^2}{4n}]/n}}{1+\\frac{z^2}{n}}$$\n",
    "\n",
    "where:\n",
    "* $\\hat{p}$ – pronounced *p-hat* – is the percentage of positive observations (opioid prescriptions in our case).\n",
    "* $z$ is the confidence level of our choosing (.95 or higher is standard)\n",
    "* and $n$ is the total number of observations (total number of prescriptions in our case).\n",
    "\n",
    "The formula looks intimidating but it's intuitive if you take the time to walk through it. Explaining the mathematics behind why this formula works is worth an entire discussion itself and is thus out of scope. So we'll just focus on applying it below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wilson_ci_lower_bound(num_pos, n, confidence=.95):\n",
    "    if n == 0:\n",
    "        return 0\n",
    "    \n",
    "    z = stats.norm.ppf(1 - (1 - confidence) / 2) # Establish a z-score corresponing to a confidence level\n",
    "    \n",
    "    p_hat = num_pos / n\n",
    "    \n",
    "    # Rewrite formula in python\n",
    "    numerator = p_hat + (z**2 / (2*n)) - (z * math.sqrt((p_hat * (1 - p_hat) + (z**2 / (4*n))) / n))\n",
    "    denominator = 1 + ((z**2) / n)\n",
    "    \n",
    "    return  numerator / denominator "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's apply this formula to our dataframe with 95% confidence, creating a new column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>num_opioid_prescriptions</th>\n",
       "      <th>num_prescriptions</th>\n",
       "      <th>opioid_prescription_ratio</th>\n",
       "      <th>wilson_opioid_prescription_ratio</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>prescriber_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>PBNDDAEOOPP</th>\n",
       "      <td>50</td>\n",
       "      <td>70</td>\n",
       "      <td>0.714286</td>\n",
       "      <td>126.310060</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PBIOORDIBSG</th>\n",
       "      <td>50</td>\n",
       "      <td>64</td>\n",
       "      <td>0.781250</td>\n",
       "      <td>114.888600</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>PDS$GG9GOGA</th>\n",
       "      <td>48</td>\n",
       "      <td>61</td>\n",
       "      <td>0.786885</td>\n",
       "      <td>109.153961</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>D4OD9B9P9O9</th>\n",
       "      <td>46</td>\n",
       "      <td>59</td>\n",
       "      <td>0.779661</td>\n",
       "      <td>105.321622</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>&amp;AGO$DEOPRR</th>\n",
       "      <td>47</td>\n",
       "      <td>58</td>\n",
       "      <td>0.810345</td>\n",
       "      <td>103.438340</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$OBDNSADD96</th>\n",
       "      <td>37</td>\n",
       "      <td>57</td>\n",
       "      <td>0.649123</td>\n",
       "      <td>101.373979</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>SEB6O6DPG3G</th>\n",
       "      <td>43</td>\n",
       "      <td>54</td>\n",
       "      <td>0.796296</td>\n",
       "      <td>95.778454</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4BOSRODA$GI</th>\n",
       "      <td>49</td>\n",
       "      <td>52</td>\n",
       "      <td>0.942308</td>\n",
       "      <td>92.095016</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>&amp;SOEB4RODOG</th>\n",
       "      <td>39</td>\n",
       "      <td>51</td>\n",
       "      <td>0.764706</td>\n",
       "      <td>90.017303</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>B$AREPGD&amp;RG</th>\n",
       "      <td>35</td>\n",
       "      <td>49</td>\n",
       "      <td>0.714286</td>\n",
       "      <td>86.151136</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               num_opioid_prescriptions  num_prescriptions  \\\n",
       "prescriber_id                                                \n",
       "PBNDDAEOOPP                          50                 70   \n",
       "PBIOORDIBSG                          50                 64   \n",
       "PDS$GG9GOGA                          48                 61   \n",
       "D4OD9B9P9O9                          46                 59   \n",
       "&AGO$DEOPRR                          47                 58   \n",
       "$OBDNSADD96                          37                 57   \n",
       "SEB6O6DPG3G                          43                 54   \n",
       "4BOSRODA$GI                          49                 52   \n",
       "&SOEB4RODOG                          39                 51   \n",
       "B$AREPGD&RG                          35                 49   \n",
       "\n",
       "               opioid_prescription_ratio  wilson_opioid_prescription_ratio  \n",
       "prescriber_id                                                               \n",
       "PBNDDAEOOPP                     0.714286                        126.310060  \n",
       "PBIOORDIBSG                     0.781250                        114.888600  \n",
       "PDS$GG9GOGA                     0.786885                        109.153961  \n",
       "D4OD9B9P9O9                     0.779661                        105.321622  \n",
       "&AGO$DEOPRR                     0.810345                        103.438340  \n",
       "$OBDNSADD96                     0.649123                        101.373979  \n",
       "SEB6O6DPG3G                     0.796296                         95.778454  \n",
       "4BOSRODA$GI                     0.942308                         92.095016  \n",
       "&SOEB4RODOG                     0.764706                         90.017303  \n",
       "B$AREPGD&RG                     0.714286                         86.151136  "
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "prescribers_df['wilson_opioid_prescription_ratio'] = prescribers_df \\\n",
    "    .apply(lambda row: wilson_ci_lower_bound(row['num_opioid_prescriptions'], row['num_prescriptions']), axis=1)\n",
    "\n",
    "prescribers_df.sort_values('wilson_opioid_prescription_ratio', ascending=False).head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Amazing! Now these are the results we would want to report. \n",
    "\n",
    "While there are plenty of prescribers with $1/1$ or $2/2$ opioid prescriptions, they now appear at the bottom of our rankings – which is intuitive. While our highest ranking prescribers have a lower `opioid_prescription_ratio` than some other prescribers, the rankings now take into account a notion of mathematical confidence.\n",
    "\n",
    "Both approaches to this problem — with or without using confidence intervals — are *technically* acceptable. However, it's easy to see that approaching the problem with a mathematical background yielded much more valuable results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# In Conclusion..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A large part of the data science job is translating open and interpretable questions into quantitative, rigorous terms. \n",
    "\n",
    "As both of these examples have demonstrated, sometimes it’s not such an easy task and many data scientists often fall short in this regard. It’s very easy to fall into the trap of using parametric correlation in a non-parametric scenario. And it’s very easy to naively sort a list of Bernoulli trials without taking into account the number of observations for each trial. In fact, it happens more often than you’d think in practice. \n",
    "\n",
    "Part of what separates a good data scientist from a great data scientist is having the mathematical background and intuition to identify and act upon situations like this. Often times what can make the difference between a solution and a great solution in data science is leveraging the right mathematical tools in the right context."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
